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SUMMARY 

Mesh deformation in response to redefined boundary geometry is a frequently encountered task in shape optimization and 
analysis of fluid- structure interaction. We propose a simple and concise method for deforming meshes defined with three-node 
triangular or four-node tetrahedral elements. The mesh deformation method is suitable for large boundary movement. The 
approach requires two consecutive linear elastic finite-element analyses of an isotropic continuum using a prescribed 
displacement at the mesh boundaries. The first analysis is performed with homogeneous elastic property and the second with 
inhomogeneous elastic property. The fully stressed design is employed with a vanishing Poisson’s ratio and a proposed form of 
equivalent strain (modified Tresca equivalent strain) to calculate, from the strain result of the first analysis, the element-specific 
Young’s modulus for the second analysis. The theoretical aspect of the proposed method, its convenient numerical 
implementation using a typical linear elastic finite-element code in conjunction with very minor extra coding for data processing, 
and results for examples of large deformation of two-dimensional meshes are presented in this paper. 

KEY WORDS: Mesh deformation, shape optimization, fluid-structure interaction, fully stressed design, finite-element analysis, 
linear elasticity, strain failure, equivalent strain, Tresca failure criterion 


INTRODUCTION 

Mesh updating in response to redefined boundary geometry is a frequently encountered task in shape 
optimization and analysis of fluid- structure interaction. The updating can be accomplished by mesh 
deformation or mesh regeneration. In contrast to a regenerated mesh, a deformed mesh retains its original 
mesh connectivity. It is advantageous to retain mesh connectivity in shape optimization and analysis of 
fluid-structure interaction. Therefore, whenever mesh deformation is applicable, it is considered a 
preferred approach. 

Gnoffo [1] and Nakahashi and Deiwert [2] employed a one-dimensional spring analogy to develop 
adaptive grid algorithms that may improve flowfleld solutions. Subsequently, several researchers [3-7] 
applied the spring analogy to obtain deformed meshes that are consistent with moving boundaries. The 
progress was from use of only linear springs to inclusion of torsional springs, from uniform spring 
stiffness to variable spring stiffness, from ignoring coupling among displacement components to the 
standard structural formulation, and from two-dimensional (2-D) applications to 3 -dimensional (3-D) 
applications. The very original spring-based method is simple and concise. However, without the aid of 
a strain tensor for multidimensional continuum, the best existing spring-based methods for 3-D 
applications do not have the advantage of simplicity or conciseness. A later advancement in mesh 
deformation techniques is the use of multidimensional linear elasticity analogy [8-12]. With the robust 
classical description of multidimensional continuum deformation, the elasticity-based methods yield more 
satisfactory results than the spring-based methods and can be formulated in a simple and concise manner. 

In some other approaches, each displacement component of a mesh movement is governed by a 
partial differential equation, e.g., Laplacian equation [13, 14] and biharmonic equation [15]. Hence, there 
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is no coupling among the displacement components. Although theses approaches are simple and concise, 
their effectiveness is limited by the independence of the displacement components. 

In the existing methods which employ variable stiffness (for the spring-based methods and elasticity- 
based methods) or variable diffusivity (for the Laplacian equation), the stiffness or diffusivity is assumed 
to be a function of measures of element (2-D or 3-D mesh cell) size (such as length, area, and volume) or 
quality (such as angle between two edges and Jacobian). For example, elements with initially smaller size 
or worse quality are assigned initially larger stiffness. However, smaller or worse elements do not 
necessarily undergo worse deformation (which is best described in terms of a strain tensor for a 
multidimensional continuum) when uniform stiffness is used. The existing methods may result in over- 
stiffening or under-stiffening of some elements, which has a negative impact on the effectiveness of mesh 
deformation. One may infer from the stiffness adjustment assumption that the existing methods tend to 
“assimilate” neighboring elements rather than tend to preserve the desirable original differences among 
the elements. 

In this paper, we present a mesh deformation method which is analogous to the Fully Stressed Design 
(FSD), an elementary design method for linear elastic structures. The objective of the present approach is 
to “assimilate” deformation of neighboring elements (not elements themselves) and preserve the original 
differences in size and quality among mesh cells. Therefore, an analogue is to design a structure of which 
“all” the elements are stressed to an allowable equivalent strain when the structure is subjected to a 
boundary movement. Obviously, the FSD is a straightforward method to get a “quick” solution to the 
pseudo design problem. It is noted that any structural design process involves a material failure model in 
addition to the linear elastic stress-strain relation and equilibrium equation that the elasticity-based 
approaches require. For the present approach, the required failure model is a strain failure model that is 
defined in terms of an equivalent strain. 

The present method is applicable to meshes defined with three-node triangular or four-node 
tetrahedral elements. The method requires two consecutive finite-element analyses of an isotropic 
continuum, using a prescribed displacement at the mesh boundaries to move the internal grid points. The 
first analysis is performed with homogeneous elastic property and the second with inhomogeneous elastic 
property. The FSD is employed to calculate, from the strain result of the first analysis, the element- 
specific Young’s modulus for the second analysis. To this end, the Poisson’s ratio of the elastic 
continuum is set to zero and a proposed form of equivalent strain (modified Tresca equivalent strain) is 
used. The formulation of the present method is simple and concise, its numerical implementation is 
convenient, and the method is effective in handling large boundary movement. Theses advantages are 
obvious in the following sections, where we present the theoretical aspect of the proposed method, its 
numerical implementation using a typical linear elastic finite-element code in conjunction with very 
minor extra coding for data processing, and results for examples of large deformation of 2-D meshes. 


THEORETICAL ASPECT 

Without loss of generality, we consider a meshed domain V as shown in Fig. 1. The boundary S 
consists of an inner boundary S t and an outer boundary S Q . The mesh is expected to deform to a new one 
that has acceptable mesh quality and meets certain boundary displacement requirements, which define 
new boundary geometry. As in the elasticity-analogy approach, the original mesh is viewed as a finite- 
element mesh for an isotropic linear elastic material subject to the boundary displacement requirements of 
the mesh deformation problem and the desired deformed mesh is defined by a finite-element solution for 
the nodal displacement of the elasticity problem. A basic assumption of the proposed method is that all 
the inhomogeneous boundary conditions for the elasticity problem are of displacement type; namely, 
there is no inhomogeneous boundary condition of traction type. 

If all the homogeneous boundary conditions are of displacement type as well, the volume average of 
the linear strain tensor depends on the displacement specified over the boundary but it is independent of 
the elastic property. This fact is obvious from the following derivation: 
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Boundary: S= S t + S o prescribed on S 

Figure 1. Schematic illustration of a linear elasticity problem with all the inhomogeneous boundary conditions of displacement 
type. 

— \e..dv = — f— (u. .+u. )dv = -^—[(un.+u.n)ds (1) 
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where i = 1, 2, or, 3 , j = 1, 2, or, 3 , u t and £ tj are the Cartesian tensors of displacement and strain, 
respectively, V represents the domain occupied by the elastic continuum, S is the boundary of F, n i is the 
outward unit vector normal to S , and F 0 is the volume of F. Because the average strain for a specific 
mesh deformation problem is fixed, a significant reduction in a strain concentration is usually 
accompanied by a strain increase in a low strain region. 

In this paper, we primarily consider 2-D meshes defined with three-node triangular elements and 3-D 
meshes defined with four-node tetrahedral elements. The mesh deformation is accomplished by two 
consecutive linear finite-element analyses using the original mesh. The first one is performed with 
uniform elastic property to produce a baseline solution followed by the second one with non-uniform 
elastic property. The nodal displacements resulting from the second analysis define the desired new 
locations of the nodes. The fully stressed design is employed to determine suitable material stiffness for 
each element in order to have more uniform deformation. The FSD is typically applied to design of a 
structure whose primary response is governed by a linear constitutive relation between a scalar strain (or 
scalar generalized strain) and a scalar stress (or scalar generalized stress). For the case of a truss, the 
scalar linear relation is £ = (1 /K)P , where £ and P are the axial strain and axial load in a rod of the truss, 
respectively, and K = AE with A and E being the cross-sectional area of the rod and Young’s modulus of 
the material of the rod, respectively. The FSD sizes the stiffness property of the aforementioned linear 
relation, so that the strain would be uniformly at its allowable value if the stress were not affected by the 
sizing. The FSD uses strain response calculated with assumed stiffness property values for an initial 
design to determine new stiffness property values for an improved design. For the case of a truss, the 
“FSD equation” used to calculate the improved stiffness of a rod is K new = {s old /£ allow )K old , where £ 0 u and 

K 0 i d are the strain and stiffness of the rod of the initial design, respectively, £ a u ow is the allowable strain, 
and K new is the stiffness estimate for the improved design. 
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To apply the FSD to 2-D or 3-D isotropic continuum under multi-axial straining without ignoring any 
tensorial strain component in the calculation, some manipulation is necessary. First, we set the Poisson’s 
ratio of the elastic continuum to zero. The setting reduces the stress-strain relation from 
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where cr. is the Cartesian stress tensor. Note that the material compliance matrix is degenerated into a 

scalar and the material with the vanishing Poisson’s ratio behaves as if it were a simple one-dimensional 
spring characterized by the same elastic constant for any type of loading. Secondly, an equivalent strain 
F(e y ) , which is a nonnegative, homogeneous, and isotropic function of the strain components, is used as 

a scalar measure of the mesh deformation. Substituting Eq. 3 into F(^) and use of the functional 
homogeneity yield 


F(s s ) = jF(a„) 


(4) 


Equation 4 is analogous to the aforementioned linear constitutive relation for rods. By the analogy, the 
“FSD equation” used to calculate the improved material stiffness for an element is 


F(e° ,d ) 


£ , = 


^ old 


(5) 


where sf 1 is the element strain result of the first analysis, E 0 i d and E new are the Young’s moduli used for 

the element for the first and second analyses, respectively, and F a u ow is the targeted limit value of the 
equivalent strain. The two types of finite elements considered in this paper are the so-called constant- 
strain elements; there is no variation in s°J d within each element. It should be noted that uniformly 

scaling up or down the Young’s modulus affects the stress solution but does not affect the strain solution 
to the elasticity problems where all the inhomogeneous boundary conditions are of displacement type. 
For convenience, the numerical value of E 0 i d (i.e., the uniform Young’s modulus for the first analysis) is 
set to unity whatever unit of stress is chosen for the analysis. 

Equation 5 cannot be employed directly to calculate E new because E new calculated with Eq. 5 for 
elements with a vanishing equivalent strain is zero and causes the singularity problem. Let F max and F min 
denote the maximum and minimum of the equivalent strain over all the elements that are not fully 
constrained by the boundary conditions, respectively. Note that 0 < F min because the equivalent strain is 
nonnegative. We only consider nontrivial problems, i.e., problems exhibiting a variation in strain; 
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therefore, F min < F max . We further restrict our investigation to problems for which the domains of analysis 
are sufficiently dominated by low-strain regions. We postulate that for these problems, F min is a 
reasonable choice for the targeted strain level, if 0 < F min and the maximum of E new calculated with the 
FSD equation is not so much larger than the minimum as to cause the singularity problem. If F min = 0 or 
(Fmax/ Fmin) E o id * s f ar t°° l ar g e to avoid the numerical problem, a numerically acceptable value should be 
chosen as the upper bound of the E new / E old ratio. It is a frequent observation that mesh deformation 
results for problems with a high F max / F min ratio is insensitive to the chosen upper bound of the E new / E old 
ratio provided that it is beyond a problem-dependent threshold. Guided by the above considerations, we 
propose the following equations to calculate the element-specific Young’s modulus for the second 
analysis: 


c = 



E nm ={\ + cT)E 0 


T = 

F -F 

max min 
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where E old = 1 with the unit of stress chosen for the analysis and 0 < c max is the only numerical parameter. 
Equations 6 and 7 imply that elements with F(s°J d ) = F min are assigned the minimum of E new , which is 
E oid , and elements with F (s ° ld ) = F max are assigned the maximum, which is either (F max / F min )E old or 

(1 + c max )E old . The numerical parameter c max should be set to as large a number as possible, subject to the 
condition that it should not render the finite-element equations prone to the singularity problem. 
Typically, c max is not less than 10 3 . It is straightforward to show that for the case of 
(Fmax - F min)/ F min^ c m a X ’ E ti- 6 is identical to the FSD equation with F allow = F min . For the case of 
c max - ( F max ~ F min) I F min or F min = 0 > it can be shown easily that the difference in E new between Eq. 6 and 
the FSD equation with F allow = (F max -F min )/c max is negligible if 0.1 < T < 1 and 10 3 < c max . Equations 6-8 
are modified FSD equations that are intended for the mesh deformation problems under present 
consideration. 

The modified Tresca equivalent strain is proposed herein, which is nonnegative, homogeneous, and 
isotropic functions of the strain components. For 3-D cases, the modified Tresca equivalent strain is 
defined by 


Fj (Sy ) — Max | v£ i £ m , c£ j | (9) 

For 2-D cases, it is defined by 

F r (s t j ) = Max{rs, - s n , es,\ (10) 

In the above, £ m < £ n < £ d are the principal strains, and 0 < r < 1 and 0 < e are numerical parameters. 
The r£ I -£ III term of Eq. 9 and the r£ I -£ II term of Eq. 10 account for failure due to combined 
maximum shear and contraction (or dilatation) on the principal planes, while the e£, term accounts for 
maximum stretching failure. The development of the modified Tresca equivalent strain and its associated 
strain failure criterion, and the determination of the two parameters are presented in the Appendix. As 
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explained in the Appendix, reasonable values such as r = 0.25 and e = 0.1 can be chosen easily though 
heuristically. The present interest is to demonstrate that the above values, which are chosen with common 
engineering judgment, are good for a number of problems rather than that certain optimum values, which 
have to be determined by a sophisticated way, may produce the best result for a specific problem. 


IMPLEMENTATION 

The present method can be implemented with a basic finite-element code, which employs constant- 
strain elements, deals with isotropic linear elastic materials, allows use of element-specific elastic 
property, and calculates nodal displacements and element principal strains. Such a finite-element code is 
virtually ubiquitous. There is only one code that is specific to the present method and needs to be 
developed, and this code is rather short. For each element, the special code uses the principal strains from 
the first analysis to calculate the equivalent strain according to either Eq. 9 or Eq. 10 and then the 
element-specific Young’s modulus according to Eqs. 6-8, and finally outputs the corresponding elastic 
material definition and its link to the element in preparation for the second analysis. If the chosen finite- 
element code allows an algebraic definition of the Young’s modulus as a material function of some 
element-dependent field variables, e.g., element temperature, then an alternative is to have the special 
code only output the value of T (Eq. 7) and its associated element identity and directly define the material 
function (Eq. 6) in the analysis model. This alternative may considerably reduce the storage of the 
analysis model. A remark about the strain result of the first analysis is worth noting. Many finite- 
element codes or post-processors generate the so-called smoothed or averaged strain. The smoothed 
strain does not exactly represent element deformation. Therefore, the smoothed strain is not a good data 
source for the calculation of the element-specific Young’s modulus. In fact, numerical findings have 
corroborated that the smoothed strain pronouncedly impairs the effectiveness of the present method. 

The implementation of the proposed method simply consists of three steps: Step 1, performing the 
first round of finite-element analysis using E old = 1 with the unit of stress chosen for the analysis; step 2, 
executing the above-mentioned special code to incorporate the element- specific E new into the analysis 
model; step 3, performing the second round of finite-element analysis using the updated model. The 
nodal displacements resulting from the second analysis define the new nodal locations, which together 
with the unchanged mesh connectivity, define the deformed mesh. 

In the present investigation, finite-element analysis is performed with the commercially available 
program NASTRAN Version 2005 of the MSC. Software Corporation. In addition, the required special 
code computes and outputs data of 2-D element temperature in a NASTRAN input format and an 
isotropic linear elastic material with the temperature-dependent Young’s modulus is readily defined 
within a few lines in the NASTRAN input file for the second round of analysis. All computations for the 
numerical examples presented in the following sections are performed on a computer with an Intel Xeon 
3.07 GHz CPU, 3 GB RAM, and the Linux operating system. 

The parameters r = 0.25, e = 0.1, and c max = 10 6 are used throughout the numerical examples 
presented in this paper. Using high-precision computing, NASTRAN should be able to handle the 
difference in element stiffness caused by c max =10 6 for any practical problem. As mentioned in the 
preceding section, we intend to evaluate the proposed method with a conveniently obtained single set of 
values of r and e for a variety of problems. 


TWO-DIMENSIONAL EXAMPLE 1 

We choose a numerical example reported in [15] as the first example in this paper. Shown in Fig. 2 is 
a 2-D mesh for a square domain whose edges are one unit long. Being similar to that used in [15], the 
mesh has 2048 elements and 1089 nodes. Both displacement components are set to zero at the lower edge 
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of the square. There is neither displacement in the x direction nor traction in the y direction at both the 
left and right edges. The x displacement at the upper edge is also set to zero. The only inhomogeneous 
boundary condition of this mesh deformation problem is a sinusoidal distribution of y displacement with 
an amplitude of 0.25 unit and a wave length of 1 unit enforced at the upper edge. The valley of the 
sinusoidal distribution is located at the midpoint of the upper edge. Each finite-element analysis for this 
example takes less than one second. 



Figure 2. Undeformed mesh of the example 1. 


The deformed mesh obtained from the first analysis is shown in Fig. 3. Although deterioration of the 
mesh is notable in Fig. 3, no grid line crossover has occurred. This suggests that the problem is still a 
mild one. However, the result obtained by solving Laplacian equations, reported in [15], indicates mesh 
penetration. 




a) b) 

Figure 3. Deformed mesh of the example 1, obtained from the first analysis: a) Entire mesh, b) Mesh around the midpoint of the 
upper edge. 


The c value computed for the second analysis is c = 3.044 xlO 2 and the mesh obtained from the 
second analysis is shown in Fig. 4. Obviously, the original mesh quality is well preserved after the 
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deformation. The result obtained by solving biharmonic equations, reported in [15], is similar to that 
shown in Fig. 4; however, the deformation is enforced in ten increments in [15]. 




a) b) 

Figure 4. Deformed mesh of the example 1, obtained from the second analysis: a) Entire mesh, b) Mesh around the midpoint of 
the upper edge. 


TWO-DIMENSIONAL EXAMPLE 2 

As shown in Fig. 5, a cylindrical object is stowed in a side bay. Three positions of the side-bay door 
are indicated in the figure. The door is slightly opened to the position A, opened to the halfway position 
B, and fully opened to the position C. The cylindrical object may be displaced to its ready position if the 
door is fully opened. The sub-domain, indicated in Fig. 5, with the cylindrical object stowed and the door 
in the position B is meshed, and the mesh is shown in Fig. 6, which has 353 elements and 211 nodes. 
Note that the thickness of the side -bay door is not modeled. Three cases of mesh deformation are 
investigated in this example: 1) mesh deformation due to door closing form the position B to the position 
A, 2) mesh deformation due to door opening from the position B to the position C; 3) mesh deformation 
due to door opening from the position B to the position C and the displacement of the cylindrical object 
from the stowage to the ready position. The external boundary is held stationary for each case. 

The coarse mesh shown in Fig. 6, which may be called a “geometry mesh,” is supposed to represent 
the boundary geometry and boundary movement with sufficient accuracy and is intended to be deformed 
by the present method, but not to be used for actual engineering analysis (e.g., fluid dynamic 
computation). Another much finer mesh, which may be called an “analysis mesh,” should be used for the 
actual engineering analysis. Once the deformation of the geometry mesh is completed, most grid points 
of the analysis mesh can be relocated by interpolating the displacement field defined on the geometry 
mesh. By using a coarse geometry mesh, both computational time and strain concentration captured in 
finite-element analysis can be reduced; therefore, applicability of the mesh deformation method can be 
increased. 

The deformed meshes obtained from the first analysis are shown in Figs. 7-9. The meshes are not 
usable because some elements around the free end of the door are crushed in all the three cases and 
several elements around the portion of the cylindrical boundary facing the door are also crushed in the 
first and third cases. 
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Figure 5. Problem definition of the example 2 



Figure 6. Undeformed mesh of the example 2: a) Entire mesh, b) Mesh around the cylindrical object. 


The respective c values computed for the second analysis for the three cases are c = 1.307 xlO 3 , 
c = 7.549 xlO 3 , and c = 1.458 xlO 2 , and the meshes obtained from the second analysis are shown in 
Figs. 10-12. It can be seen in the figures that the inhomogeneous elastic property computed by the 
proposed method conduces to far more uniform deformation and thus substantially improves the mesh 
quality. For example, comparing Fig. 9a with 12a Fig, one may see much more notable propagation of 
deformation toward the low strain regions due to the introduction of the property variation. No grid line 
crossover occurs in any of the deformed meshes obtained from the second analyses. Note that each of the 
three deformed meshes is derived from the same original mesh (Fig. 6) by using one increment. 
Therefore, a mesh for any intermediate instance during the process of changing the geometry, from the 
case 1 to the case 2 and through the case 3, can be conveniently derived from the same original mesh 
using one increment as well. 
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a) b) 

Figure 7. Deformed mesh, door closing, the first analysis: a) Entire mesh, b) Mesh around the cylindrical object. 




Figure 8. Deformed mesh, door opening, the first analysis. 




a) b) 

Figure 9. Deformed mesh, door opening, the object displaced to the ready position, the first analysis: a) Entire mesh, b) Mesh 
around the cylindrical object. 
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a) b) 

Figure 10. Deformed mesh, door closing, the second analysis: a) Entire mesh, b) Mesh around the cylindrical object. 



Figure 11. Deformed mesh, door opening, the second analysis. 



Figure 12. Deformed mesh, door opening, the object displaced to the ready position, the second analysis: a) Entire mesh, b) 
Mesh around the cylindrical object. 
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TWO-DIMENSIONAL EXAMPLE 3 


Shown in Fig. 13 is a 2-D unstructured mesh intended for a transonic flow over an airfoil. The mesh 
has a total of 18812 elements and 9598 nodes with many of them clustering near the airfoil surface. The 
chord of the airfoil is parallel with the x-axis and has a length of 1 unit. To demonstrate the effectiveness 
of the present method, we attempt to deform the mesh in four different ways: 1) by a translation of the 
airfoil by 0.5 unit in the flow (x) direction, 2) by a 0.5 -unit upward (y) translation, 3) by a 45° rotation 
about the mid-chord point, and 4) by cambering the airfoil shape via imposing an arc-shaped distribution 
of upward (y) displacement (zero at both leading and trailing edges, and a peak displacement of 0.125 unit 
at the mid-chord point). The camber movement is often encountered in a design optimization process, 
while the rigid-body motion cases may appear in both design optimization and unsteady aerodynamic 
analysis. Each NASTRAN run takes about five seconds. 

Figures 14-17 show the mesh patterns resulting from the first round of finite-element analysis. The 
uniform elastic property is employed at this stage and the resulting meshes are of very poor quality. It 
appears that the near-wall mesh points have been dragged along with the airfoil motion. In particular, 
grid line crossover occurs at the trailing edge for all the cases. Furthermore, the near- wall elements lose 
the quasi-orthogonality of the baseline mesh. 

The respective c values computed for the second analysis for the four cases are c = 6.922 xlO 2 , 
c = 7.009 x 10 2 , c = 4. 156 x 10 3 , and c = 2.660 x 10 4 , and the meshes obtained from the second analysis 
are shown in Figs. 18-21. Evidently, the artificial temperature fields created for the four test cases and the 
temperature-dependent elastic material have been effective in handling large boundary movement. The 
leading and trailing edges are problem-free and no grid line crossover is observed in the results. For the 
two translation cases, the orthogonality of the near- wall mesh is almost completely preserved. For the 
other two cases, the near-wall mesh, although less orthogonal than for the translation cases, still has good 
quality. The variable elastic property plays a crucial role in achieving the desirable mesh deformation. 
The impact of the property variation can be visualized, for example, by comparing Fig. 16a (invalid 
mesh) with Fig. 20a (valid mesh). It is seen from the comparison that the deformation is successfully 
propagated toward the low strain regions. A consequence of the propagation is the dramatic reduction in 
strain concentration, as indicated by Figs. 16b and 20b. 


CONCLUDING REMARKS 

We extend the idea of elasticity analogy to the idea of FSD analogy. We set the Poisson’s ratio to 
zero in order to retain all the strain components in the FSD formulation and thus obtain the formulation 
suitable for general mesh deformation problems. The strain failure model based on the modified Tresca 
equivalent strain, which we propose specifically to characterize mesh deformation, is used as the failure 
criterion required by the FSD. The entire formulation of the method is concise and can be conveniently 
implemented by using a typical linear elastic finite element code. The application of the method to the 
various large deformations of the 2-D meshes is demonstrated, where consistently satisfactory mesh 
quality is obtained by using the conveniently obtained single set of values of the numerical parameters, 
i.e., r = 0.25 and e = 0.1. In contrast, the existing methods from the literature exhibit difficulties in 
solving similar problems. 

The present method may be easily adapted to 2-D meshes containing four-node quadrilateral elements 
and 3-D meshes containing eight-node hexahedral elements. Because generally, strain distribution in 
such elements is not constant, the equivalent strain needs to be evaluated by using the actual strain at each 
corner of each element and then the maximum equivalent strain for the element is determined. The 
maximum equivalent strain for each element should be used as F(s °! d ) for the element in the above 
formulation of the proposed mesh deformation method; otherwise, the formulation remains the same. 


12 




Figure 13. Undeformed mesh of the example 3: a) Entire mesh, b) Mesh around the airfoil, c) Mesh around the leading edge, 
d) Mesh around the mid chord, e) Mesh around the trailing edge. 


13 




Figure 14. Mesh deformed by a 0.5-unit translation of the airfoil in the flow (x) direction, obtained from the first analysis: a) 
Entire mesh, b) Mesh around the airfoil, c) Mesh around the leading edge, d) Mesh around the mid chord, e) Mesh around the 
trailing edge. 
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Figure 15. Mesh deformed by a 0.5-unit upward (y ) translation of the airfoil, obtained from the first analysis: a) Entire mesh, b) 
Mesh around the airfoil, c) Mesh around the leading edge, d) Mesh around the mid chord, e) Mesh around the trailing edge. 
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Figure 16. Mesh deformed by a 45° rotation of the airfoil about the mid-chord point, obtained from the first analysis: a) Entire 
mesh, b) Mesh around the airfoil, c) Mesh around the leading edge, d) Mesh around the mid chord, e) Mesh around the trailing 
edge. 
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Figure 17. Mesh deformed by cambering the airfoil shape, obtained from the first analysis: a) Entire mesh, b) Mesh around the 
airfoil, c) Mesh around the leading edge, d) Mesh around the mid chord, e) Mesh around the trailing edge. 
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Figure 18. Mesh deformed by 0.5-unit translation of the airfoil in the flow (x) direction, obtained from the second analysis: a) 
Entire mesh, b) Mesh around the airfoil, c) Mesh around the leading edge, d) Mesh around the mid chord, e) Mesh around the 
trailing edge. 
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Figure 19. Mesh deformed by a 0.5-unit upward (y ) translation of the airfoil, obtained from the second analysis: a) Entire mesh, 
b) Mesh around the airfoil, c) Mesh around the leading edge, d) Mesh around the mid chord, e) Mesh around the trailing edge. 
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Figure 20. Mesh deformed by 45° rotation of the airfoil about the mid-chord point, obtained from the second analysis: a) Entire 
mesh, b) Mesh around the airfoil, c) Mesh around the leading edge, d) Mesh around the mid chord, e) Mesh around the trailing 
edge. 
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Figure 21. Mesh deformed by cambering the airfoil shape, obtained from the second analysis: a) Entire mesh, b) Mesh around 
the airfoil, c) Mesh around the leading edge, d) Mesh around the mid chord, e) Mesh around the trailing edge. 
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Numerical efficiency is critical if the present method is applied to large 3-D meshes, in particular, 3- 
D meshes for viscous fluid dynamic analysis. A way to shorten computational time is to deform a coarse 
geometry mesh first, then map “off-boundary nodes” of a fine analysis mesh to the deformed state, and 
finally perform deformation of a boundary layer of the analysis mesh to displace “near-boundary nodes” 
(i.e., nodes inside the boundary layer). The near-boundary nodes are those that may not be accurately 
relocated by the mapping. To ensure that the number of the near-boundary nodes is acceptably small, 
boundary geometry and boundary movement need to be modeled with the geometry mesh sufficiently 
accurately. To achieve the accuracy and keep the number of degrees of freedom of the geometry mesh 
small at the same time, degenerated higher-ordered elements of which all the non-corner nodes are placed 
on the boundary of the geometry mesh may be used. Besides the time saving, weaker strain concentration 
is captured in finite-element analysis due to the coarseness of the geometry mesh. Results of 
implementation of 3-D geometry meshes will be reported in a sequel. 
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APPENDIX 

In the existent methods of elasticity analogy or spring analogy, mesh deformation is calculated by 
structural analysis. All the existent methods are based on the equation of equilibrium and certain stress- 
strain relations. However, none of the existent method has employed the concept of material failure, 
which is an essential element of classical structural analysis. 

A strain failure criterion can be set up to determine whether a mesh cell is so severely deformed that it 
cannot be used for computation. Typically, a strain failure criterion is a nonnegative and homogeneous 
function of the strain tensor. Such a failure criterion can be cast into a simple statement that a material 
fails when an equivalent strain reaches an allowable value. 

We borrow the notion of the Tresca equivalent stress to define a Tresca equivalent strain as the 
following: 


Ft ( £ ij ) = Tfaxj £j - £ n 


■-*/!} 


(Al) 
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where £ 7 , £ /n and s m are the principal strains. The strain failure criterion associated with F r ° is “A 
material fails if F r ° > £ f c , where 0 < 4 < 100 % is a material parameter, magnitude of the uniaxial 
compressive failure strain.” Obviously, F r ° is a nonnegative, homogeneous, and isotropic function of the 
strain tensor. The failure surface represented by = s f c is an open-ended hexagonal cylinder with 

the centerline defined by £, = s u = s m in the principal strain space. Note that the Tresca equivalent strain 
is a direct measure of the maximum shear strain, and it is independent of the hydrostatic strain 
component. These are desirable features because it is shear deformation of a mesh cell (not volumetric 
deformation) that accounts for primary deterioration of the mesh cell. However, deterioration of mesh 
quality is not exactly independent of the volumetric deformation. For example, an initially perfect mesh 
cell (say, an equilateral triangle for the 2-D case) may expand or contract to an unacceptable size without 
any change in the aspect ratio or twist. 

For brevity but without loosing generality, we assume £ m < s n < £7 in the following derivation. 
Therefore, the above equivalent strain reduces to 

Fj (£ \j ) — Max{Sj — £ n , £ n — £ [U , £j — £ In } — £ f — £ [[I (A2) 

To account for the volumetric deformation, the following modified Tresca equivalent strain is introduced: 

Fj (£jj ) — Max{a(£j — £ n ) — b{£ f + £ u ), a(£ u — £ in ) — bi y £ u + £ in ), Cl{£j — £ jjj ) — b(y£ j + £ jjj ), C£ j } (A3) 

where a , b , and e are material parameters. The strain failure criterion associated with F T is “A material 
fails if F r > £ j c .” Among the four competing entities in Eq. A3, the leading three, e.g., 
a(£ T -£ II )-b(£ I +£ n ) , account for failure due to combined maximum shear and contraction (or 
dilatation) on the principal planes, while the last one, e £ } , accounts for maximum stretching failure. We 
impose 0 < a so that the maximum shear strain on a principal plane has a positive contribution to the 
associated combined failure mode. The inequality 0 < b is imposed so as to ensure that a pure and severe 
enough contraction on any principal plane may cause failure. We also impose 0 < e so that the maximum 
stretching strain (<57 with a positive value) may cause failure. Note that without the e£j term in Eq. A3, 

F t could be negative and the failure surface defined by F T (£ ij ) = £ f c would not be a closed surface. The 
condition that the uniaxial compressive strain state (< 57 , £ in £ U1 ) =( 0 , 0 , ) is on the failure surface 

leads to a + b = l. By using a + b = 1 and £ m < £,, < £7 , we may simplify Eq. A3 to the following: 

F t (£ jj ) = Max {rs 7 - £ m , e£ l } ( A4) 

where r = a-b. With the additional parameters, r and e , we may impose two conditions on the failure 
surface. Firstly, the pure contraction strain state (£ 7 , £ in £ ni ) = (-100%, -100%, -100%) is on the failure 
surface. Therefore, 


Max{-r + 1, - ej = -r + 1 = £ f c (A5) 

from which we may infer 0 < r < 1 . Secondly, the pure expansion strain state 

(£ 7 , £ n , £ m )= (£*/, £■/, £ J ,) is also on the failure surface, where 0<£ > / is the maximum stretching 
failure strain. Therefore, 
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Max{(r- !)<?/, es{} = es{ = £ f c 


(A6) 


Once values of and £•/ are specified, satisfying o <s f c< 100% and 0 < s{ , one can easily calculate r 
and e by using Eqs. A5 and A6. One may infer from Eqs. A4 and A6 that a material cannot survive a 
maximum principal strain greater than or equal to s{ . An alternative approach is to use the uniaxial 
tensile failure strain, s f T , instead of the magnitude of the uniaxial compressive failure strain, ^ . Note 
that s f T and s{ satisfy 0 <8 f T <s{ . To implement the alternative approach, we postulate that the 
uniaxial tensile failure is a combined failure mode, namely, 

Max{r££ , e£j } = r£ f T = £ f c (A7) 

Substituting Eq. A7 into Eq. A5, one may derive the following: 

r = 1/(4 +1 ) ( A8 ) 

Substituting Eq. A7 into Eqs. A6, one may derive the following: 

e = r£y / £■[ (A9) 

Therefore, once values of £ J r and 6*/ are specified, satisfying 0<£^ < £ J , , one can easily calculate r and 
e by using Eqs. A8 and A9. The failure surface and material parameters are illustrated by Fig. Al. The 
figure helps one determine r and e in a straightforward manner. 



(-100%, -100%) 

Figure Al. Failure surface, based on the modified Tresca equivalent strain, and the associated material parameters. 

The postulate that the uniaxial tensile failure is a combined failure mode reflects a general 
observation that unsuitability of a mesh cell caused by uniaxial staining is an effect of the distortion or 
aspect ratio in addition to the cell size. On this postulate, determination of the parameter r substantially 
depends on how one perceives equivalence of the uniaxial tensile failure strain state, i.e., (£ n £ in £ ni ) = 
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(<£•/, 0, 0), and the pure contraction failure strain state, i.e., (e n s in s m ) = (-100%, -100%, -100%). For 
example, if one considers (s n s in s m ) = (300%, 0, 0) as undesirable a state of strain failure as 
(s n s in £„,) = (-100%, -100%, -100%), then r, which can be calculated with Eq. A8 and is indicated as a 
slope in Fig. Al, is equal to 0.25. With the decided values of s? and r, one may proceed to select a value 
of s{ in order to calculate the parameter e with Eq. A9. For example, if s f T = 300% is chosen, one may 
choose s{ = 750% . The selected values result in e = 0.1 . 


25 



